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A DIGITAL COMPUTER PROGRAM FOR CALCULATING THE 
PERFORMANCE OF SINGLE- OR MULTIPLE -DIAPHRAGM 
SHOCK TUBES FOR ARBITRARY EQUILIBRIUM 
REAL GAS MIXTURES 

By William L. Grose, John E. Nealy, and Jane T. Kemper 
Langley Research Center 

SUMMARY 

A computer program written in FORTRAN IV language is presented which deter- 
mines the performance of a shock tube for arbitrary equilibrium real gas mixtures. For 
specified initial gas composition and charging conditions, the program output includes 
velocity, pressure, density, enthalpy, temperature, sound speed, and species mole frac- 
tions at any point in the shock tube cycle. 

Several representative calculations illustrate the utility of the program. The pro- 
gram is applicable to both simple and buffered shock tubes as well as to the expansion 
tube. 


INTRODUCTION 

The shock tube has been widely utilized as an experimental device for producing 
high enthalpy flows for research in such areas as aerodynamic testing, chemical kinetics, 
and radiation gas dynamics. Much of the utility of the shock tube is due to the capability 
of generating a broad range of test conditions in various gas mixtures. 

Numerous current experimental investigations are being conducted in gas mixtures 
simulating planetary atmospheres that differ markedly from air. Before initiating a test 
program, it is essential to ascertain the theoretical performance of the shock tube for the 
particular gas mixture being used in order to aid in determining test conditions and in 
analyzing data. Previous shock tube performance investigations are too numerous to 
describe herein, but they generally can be divided into two classes. Both simple and 
buffered shock tubes are discussed in reference 1 and this study is an example of the 
analyses that assume a perfect gas throughout the operating cycle. In reference 2 an 
investigation of the expansion tunnel is presented and is illustrative of those analyses that 
treat the shocked gas in the driven chamber as a real gas in thermodynamic and chemical 
equilibrium, but the driver gas in the unsteady expansion process is assumed to be 



perfect. The rather exhaustive treatment of the perfect gas performance of the simple 
shock tube in reference 3 does, however, include a discussion of real gas effects. 

It is desirable to construct a digital computer program to calculate shock tube per- 
formance for arbitrary driver and driven gases so that the gas at any point in the operating 
cycle may be treated either as a perfect gas or as an equilibrium real gas. This report 
describes such a program based upon the assumptions of inviscid one -dimensional flow, 
no mixing at the interfaces, and instantaneous diaphragm bursts. The program is appli- 
cable to simple or buffered shock tubes and expansion tubes. Performance parameters 
for some typical gas mixtures and operating conditions are presented to demonstrate the 
utility of the program. 


SYMBOLS 


a sound speed, cm/sec 

h enthalpy, ergs/g 

p pressure, dynes/cm^ 

p 0 standard atmospheric pressure, 1.01325 x 10® dynes/cm2 

R universal gas constant, 8.31469 x 10^ ergs/mole-°K 

s entropy, ergs/g-°K 

T temperature, °K 

t time, sec 

u flow speed, cm/sec 

Ug shock speed, cm/sec 

x distance, cm 

y ratio of specific heats 

ju molecular weight 
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p 


density, g/cm3 


Subscripts: 

0 reference conditions 

1 state of quiescent gas in front of normal shock 

2 state of gas behind normal shock 

3 state of expanded driver gas 

4 initial driver gas conditions 

5 state of test gas in expansion tube 

10 state of quiescent gas in double -diaphragm shock tube 

ANALYSIS 

A distance-time diagram for the simple shock tube is schematically depicted in fig- 
ure 1(a). Figure 1(b) illustrates the quiescent driver and test gases prior to diaphragm 
burst at t = 0. Figure 1(c) indicates the wave system at t = t a . The test gas in 
state (T) is compressed and heated by the normal shock wave to state (2), while the 
driver gas in state (4) undergoes an isentropic unsteady expansion to state (IT). 

In practice the shock tube is often designed to operate at a fixed driver pressure 
and the driven chamber pressure (and the buffer chamber pressure for the buffered shock 
tube) is adjusted to produce the desired test conditions. For specified initial charging 
conditions (T) and (4), the unsteady expansion must be solved for state ( 3 ) and the nor- 
mal shock wave must be solved for state . The respective solutions for states (2) 
and ( 3 ) are compatible only when the pressures and velocities of the two states are equal 
(i.e., p 2 = P3 and u 2 = u 3 ). 


Unsteady Expansion 

The governing differential equation (ref. 4) for an isentropic unsteady expansion is 


du = 



( 1 ) 
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where the negative and positive signs refer to upstream and downstream waves, respec- 
tively. Then, the velocity increment imparted to the flow by the expansion is 


Au = Ug - U4 = + 



( 2 ) 


In order to relate the pressure and velocity in the as yet undetermined state (3) for 
an equilibrium real gas, it is necessary to resort to numerical techniques to evaluate the 
integral in equation (2) . The analysis formulated herein utilizes the equilibrium program 
of reference 5. The program includes ionization and dissociation and involves the fol- 
lowing assumptions: 

(1) The mixture is composed of ideal gases. 

(2) For diatomic species the rigid-rotor harmonic oscillator model is used with 
vibrational- rotational corrections for each electronic configuration. 

(3) Only electronic levels with principal quantum number less than or equal to 5 are 
included. 

The free energy of each of the species is determined from the partition function of quan- 
tum statistical mechanics. The equilibrium composition is then arrived at by minimiza- 
tion of free energy. 

This program is used to generate an array of the thermodynamic state variables 
p, h, a, and s for the gas mixture over selected ranges of pressures and tempera- 
tures. Interpolation within the array permits tabulation of p, p, h, T, and a for 
constant s = S4. The integral in equation (2) can then be evaluated by a second-order 
Gaussian quadrature (ref. 6) between limit I14 and an assumed limit 113. The pres- 
sure P3 is then determined from the tabulation at the assumed I13. It is thus possible 
to determine a unique correspondence between pg and U3 at constant s^. 

For a perfect gas, equation (2) can be integrated in closed form to yield 



Normal Shock Wave 

The conservation equations valid across a moving normal shock wave are 


p l U S = p 2( U S - u 2) 
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(5) 


Pi + PiU s 2 


P2 + ^2 ( U S " u 2)' 


hi 




( 6 ) 


The solution to equations (4) to (6) for an equilibrium real gas requires an iterative 
method to determine the pressure and velocity in state (2). This analysis utilizes the 
normal shock program of reference 7, which is based upon a free-energy minimization 
technique for calculating equilibrium thermodynamic state variables coupled with the con- 
servation equations and a modified Newton-Raphson iteration method. The equilibrium 
thermochemical calculations involve the same assumptions as those in reference 5. For 
a specified state (T) and incident shock velocity, this method determines P 2 and U 2 < 
For a range of values of Ug, the correspondence between P 2 and U 2 can then be 
established. 

For a perfect gas, the following relation between P 2 and U 2 can be readily 
derived: 


u 2 = 



r \Pi 


_2z_ 

y + 1 

yPl y + 1 ) 


1/2 


(7) 


Method of Solution 

The solution generated by the program described herein is based upon the previ- 
ously noted requirement that pg equal P 2 and ug equal ug. By following the method 
outlined for the solution of the unsteady expansion, a monotonically decreasing sequence of 
pressure as a function of velocity can be calculated. In similar fashion, the normal shock 
solution yields a monotonically increasing sequence of pressure as a function of velocity. 
The method of false position coupled with a Lagrangian interpolation (ref. 6) is then used 
to determine the crossover point of the two sequences. The crossover values of p and 
u are the required solution. 


PROGRAM DESCRIPTION 

For application to the simple shock tube (fig. 1), the program contains the four 
options corresponding to the combinations of either a real or perfect driver gas and 
either real or perfect driven gas. For the perfect gas option, it is only necessary to 
solve the closed-form equations (3) and (7) for the expansion and normal shock, 
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respectively. For the perfect gas expansion, the required input consists of the initial 
conditions p 4 , T 4 , and u 4 and the constants p 4 , y 4 , and R; in addition, a range of 
velocities 


u 4 


< 


u 3 < u 4 + 


2 \l yRT 4 

r- iy m 


is required to generate a velocity-pressure ( U 3~P 3 ) curve. 

For the perfect gas normal shock, the required inputs are pj, Tj, ju^, R, yj, 
and a range of p 2 values 


Pi < P 2 < P 4 


from which a pressure-velocity (P2"U 3 ) function is computed. 

The real gas expansion calculation utilizes the program of reference 5 to generate 
thermodynamic data. In addition to the input required for the program of reference 5, 
suitable ranges of temperatures and pressures must be included to define the array of 
thermodynamic variables. For the constant entropy value s 4 , the program reads from 
the array h as a function of 1/a and then performs the numerical integration 


u 3 


- u 4 = 



s 4 


from which the pressure- velocity (p 3 -u 3 ) function is determined. The option exists for 
readout of velocity-pressure (u 3 -p 3 ) curves for various entropy values, as well as U3 
as a function of h and T, or for proceeding directly to the solution of the normal shock. 


The program of reference 7 is used to determine equilibrium conditions behind a 
normal shock for a range of Ug values and for given p 4 , Ti, and gas composition. 
Each Ug value thus leads to a value of p 3 and corresponding value of u 3 , which 
define the pressure-velocity function to be matched to the velocity-pressure (u 3 -p 3 ^ 
curve from the expansion calculation. 


For expansion tube or buffered shock tube performance calculation, the program is 
applied twice to the simple shock tube cycle (fig. 2). The conditions behind the normal 
shock (state (2)) in the first calculation become the driver gas conditions (state (T)j for 
the second calculation. It should be noted that the driver velocity for the second calcula- 
tion is no longer zero, since instantaneous diaphragm rupture has been assumed (i.e., no 
reflection). The fact that the gas behind the first incident shock in such a system may be 
dissociated and ionized emphasizes the necessity for a real gas calculation of the expan- 
sion phase. 
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The program FORTRAN listing and a flow chart are presented for reference in the 
appendix. 


SAMPLE CALCULATIONS 


Several sample calculations illustrate the application of the program to the deter- 
mination of performance for single- and double -diaphragm shock tubes. The driver pres- 
sures chosen are such that intermolecular forces may be neglected and the assumption of 
a mixture of ideal gas is valid. 

Figure 3 shows the velocity-pressure (U3-P3) curve for the equilibrium real gas 
isentropic expansion of helium and the velocity-pressure (U2-P2) curves for normal shock 
waves in equilibrium air for four initial driven chamber pressures with Tj = 300° K. 

The helium driver conditions chosen are T4 = 15 000° K and P4^P 0 = 315. The cross- 
over point of the curves determines P3 = P2 and U3 = U2 for each driver-driven con- 
figuration. From these results the shock speeds corresponding to each calculation may 
be evaluated. The performance plot of Ugj as a function of Pijp 0 for these driver 

conditions is given in figure 4. 


Figures 5 to 9 illustrate a typical set of performance calculations for the expansion 
tube. This type of facility (ref. 8) consists of three sections separated by two destructible 
diaphragms as shown in figure 2. One section contains the driver gas ((4))> the interme- 
diate section contains the test gas and the remaining section contains the acceler- 

ating gas ((l0)) • 

For the performance calculations, the simple shock tube cycle was applied first to 
regions (4) and (7) in order to obtain conditions in region (2) . Then the calculation 
was repeated, with conditions in region (2) as the driver and the gas in region (H^) as 
the driven gas, which led to the determination of the test gas conditions in region ^5^). 

The following initial conditions were chosen for these computations: 


= 30; T4 = 2500° K; helium (perfect gas) 

10~ 3 g 10" 1 ; T x = 300° K; air (real gas) 
p o 

10'3 < K)- 2 ; T 10 = 300° K; helium (perfect gas) 

Pn 


In this calculation the temperatures experienced by the helium were not great enough to 
cause a noticeable departure from perfect gas behavior. Figure 5 shows the shock speed 
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in air Ug ? i as a function of Pj^Pq resulting from the first diaphragm burst. The final 

shock speed in helium in the accelerating chamber is shown in figure 6 for = 5 X 1CT 4 , 

Po 

10“^, and 5 x 1CT 2 . The conditions of most interest for this facility are those of the gas 
in region (5), which has been shocked in the first cycle and expanded and cooled in the 
second cycle. The test gas conditions which may be obtained from the given initial con- 
ditions are shown in figures 7, 8, and 9. As would be expected, the families of curves for 
shock speed Ug iq anc * test gas speed ug are similar (figs. 6 and 7). The performance 
data indicate that the facility may be operated in a manner in which the test gas speed and 
pressure do not change appreciably although the temperature may vary from several 
thousand degrees to several hundred degrees Kelvin. 

This program was also used to make a set of calculations for the 3.8-inch double- 
diaphragm shock tube at the Langley Research Center in order to predict the real gas 
performance for a test gas mixture of 90 percent N2 and 10 percent CO2. The calcula- 
tions for a double -diaphragm, or buffered, shock tube are carried out in the same manner 
as those described for the expansion tube. However, in this facility the test gas is con- 
tained in the third chamber instead of in the intermediate section; hence, the results of 
primary interest are the shock speeds in the test gas for a range of Pi^/P 0 and 
which are shown in figures 10 and 11. 

Figure 10 gives the performance of this shock tube for the following initial 
conditions: 




T4 = 300° K; helium (perfect gas) 


5 x 10" 2 



Tj = 300° K; helium (perfect gas) 


5 X 10 -4 


<£i2 < 5 x 10" 2 ; 
Po 


Tjq = 300° K; 90 percent N2 and 10 percent CO2 (real gas) 


Figure 11 shows the performance for these same initial conditions for gases in 
regions (T) and (To) but with perfect hydrogen used instead of helium for gas in 
region (T). It has been observed experimentally that for the initial conditions for which 
this facility is operated, the performance calculations give a reasonably accurate pre- 
diction of actual behavior. 

The results of all sample calculations were checked where possible with the results 
of references 1, 2, 4, 7, and 8. In all calculations agreement was within 1 percent. 
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resume: 


A computer program written in FORTRAN IV language is presented which deter- 
mines the performance of a shock tube for arbitrary equilibrium real gas mixtures. For 
specified initial gas composition and charging conditions, the program output includes 
velocity, pressure, density, enthalpy, temperature, sound speed, and species mole frac- 
tions at any point in the shock tube cycle. 

Several representative calculations illustrate the utility of the program. It was 
noted that these calculations agreed with results of other investigations. 

It may be concluded that the program has both generality and flexibility for com- 
puting shock tube performance in arbitrary gas mixtures. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., June 10, 1968, 

129-01-03-10-23. 
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APPENDIX 


SHOCK TUBE PERFORMANCE PROGRAM LISTING AND FLOW CHART 

The FORTRAN listing and a schematic flow diagram of the program are presented. 
The program is composed of nine subroutines as follows: 

(1) PEREX — Computes pg as a function of Ug for a given p^, U4, T4, and R 

(2) ROGO — Computes an array of thermodynamic state variables for a given range 

of p and T (program of ref. 5) 

(3) INTER - Uses a three-point Aitkens interpolation formula to find p, h/RT, 

and a/a 0 at constant (s/R )4 

(4) INTEG — Computes ^ ~ by Gaussian quadrature 

(5) PERNS — Computes U2 for a given pj, P2, Tj, yj, and R 

(6) NORMAL — With given p, Tj, and u, computes normal shock properties of 

gas 

This subroutine makes use of an equilibrium program which computes thermo- 
dynamic properties of a gas at a given T and p (program of ref. 7). 

(7) SOLUT — Given (P2> u 2) anc * (P3> u 3)> finds solution to curves 

(8) FTLUP — Uses Aitkens interpolation formula for three points; finds Ug to 

correspond to U2 

(9) SLITE — Sets a flag to be checked later on in subroutine 

Subroutines (1) to (4) apply to the expansion phase only and the program may be 
stopped at this point. Also, ug, pg, Tg, (h/RT)g, and (a/ao)g may be punched out 
to save for later computer runs so that the curves do not have to be recomputed. Sub- 
routines (5) to (9) apply to the normal shock portion of the program. 

The program listing and flow chart are reproduced in the following pages. 
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PROGRAM D 1 282 ( I NPUT , OUTPUT % PUNCH * T APE5= I NPUT , T APE6=0UTPUT ♦ TAPE9 ) 
REAL MU 

INTEGER PUN, COMPUT. REAL .EXPO 
DIMENSION AM (30 ) ,HM (30 ) 

INPUT FOR EXPANSION - REAL GAS 

COMMON PN ( 30 ) , TM ( 30 ) , NUMT . NCAPX , S0R4 , U4 

INPUT FOR EXPANSION - PERFECT GAS 

COMMON P4 ♦ T4 , U3 ( 30 ) 

INPUT CODES 

COMMON PUN, EXPO. REAL « COMPUT 
INPUT FOR NORMAL SHOCK 

COMMON PI . T 1 »P2 (30 ) , NUMP . USTO ( 30 ) ,U2(30> 

COMMON I SPEC (30 > . UMOL ( 1 0 ) , M , N ♦ BETA ( 5 ) .NS* 

1 AMC.NB.NBTA (5) 

INPUT FOR PERFECT EXPANSION AND SHOCK 

COMMON MU. GAMMA »R 

GAS NAME 
C 

COMMON NAME ( 10 ) 

C 

COMMON/BLOCK/HM . AM , MM »P3 ( 30 ) 

C 

NAMEL I ST/EXP/PN . TM , NUMT .NCAPX . S0R4 . U4 . P4 . T4 . U3 , MU , GAMMA *R , PUN, 

1 EXPO . REAL . COMPUT 

NAMEL I ST/NORMS/PI ♦ T 1 , P2 .NUMP , USTO ♦ U2 , I SPEC . JMOL . M .N . BETA ,NS .NB , 

1 NBTA ,MU • GAMMA , R .REAL , COMPUT 
EXPO=l 

CALL SLITE(1 ) 

NS = 0 
REAL =0 
COMPUT =0 

1 READ(5. 2010) NAME 
2010 FORMAT ( 1 0 A6 ) 

READ ( 5 « EXP ) 
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M WRITE(6.2008 ) 

2008 FORMAT (31 HI SHOCK TUBE PERFORMANCE PROGRAM// 1 4H ENG* J* NEALY* 10X* 
1 1 5HPR0G • J. KEMPER) 

C 

C I F ( COMPUT =1* OR 0 READ IN U3.P3 « COMPUTE U3« P3 

C 

IF (COMPUT *EQ*0 ) GO TO 10 
C 

C READ IN U3«P3 
C 

IF(REAL.EQ*1 ) GO TO 5 
READ (5* 2005) SOR4.U4 

2005 FORMAT ( 1 4X * El 5* 8 ♦ 5X « E 1 5 *8 ) 

GO TO 6 

5 READ ( 5 * 2006 ) P4.T4 

2006 FORMAT ( 1 6X *E15*8«5X«E1S*8) 

6 READ(5.2007) ( P3 ( I ) « I = 1 «NUMT ) 

2007 FORMAT (4E17.8) 

READ (5*2007) (U3 ( I ) * I = 1 «NUMT ) 

GO TO 101 
10 CONTINUE 
C 

C IF REAL *0*1 REAL* OR PERFECT GAS 

C 

IF (REAL*EQ.O ) GO TO 15 

CALL PEREX ( P4 « U4 * GAMMA * R ♦ T4 * MU « U3 * NUMT « P3 ) 

GO TO 101 

15 CALL R0G0(PN.TM.NUMT*NCAPX.GAMMA«S0R4 ) 

MM=NUMT 

CALL 1NTEG (U3*U4 ) 

WRITE (6 *201 ) S0R4. ( 1 * TM ( I ) « P3 ( I ) ♦ HM ( 1 ) • AM ( I ) • U3 ( I ) * 1 = 1 «NUMT ) 

201 FORMAT ( 7H1 OUTPUT/ 1 5X , 4HS/R=E 1 4 • 7/3H M . 2X ♦ BX « 1 HT *8X ♦ 1 OX ♦ 1 HP « 

1 1 8X * 1 HH « 1 8X * 1 HA « 

1 18X.2HU3/ ( I 3 ♦ 5 ( 2X ♦ E 1 7 • 8 ) ) ) 

101 IF (PUN*EQ* 1) GO TO 110 

100 WRITE (6*2009 )P4 . T4 , ( U3 ( I ) * P3 < I > . I = 1 .NUMT ) 

2009 FORMAT ( 1H1 . 2X . 3HP4 = E 1 5 . 7 *2X . 3HT4 = E 1 5 .7/8X . 2HU3 * 1 7X . 2HP3/ 

1 (E17.8.2X.E17.8) ) 

102 IF (EXPO.EQ* 1 ) GO TO 150 
GO TO 1 

C 

C PUNCH U3 * P3 

C 

no if(real*eq*o> go to hi 

PUNCH 2000 *P4. T4 
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2000 FORMAT (11 HPERFECT GAS * 2X . 3HP4=E 1 5. 8 .2X . 3HT4=E1 5.8 ) 

GO TO 112 

111 PUNCH 2001. S0R4.U4 

2001 FORMAT ( 8HREALGAS . 2X » 4HS/R=E 15.8.2X. 3HU4 = E 15.8) 

112 PUNCH 2003. ( P3 ( I ) . I =1 . NUMT > 

2003 FORMAT (4E17.S.4X. 2HP3 ) 

PUNCH 2004 . ( U3 ( 1 ) , I = 1 , NUMT ) 

2004 FORMAT (4E17.8.4X.2HU3) 

C 

GO TO 102 

NORMAL SHOCK 

150 READ (5. NORMS ) 

CALL SL I TE ( 2 ) 

1515 WR 1TE(6.206) NAME .P4.T4.T1 

206 FORMAT ( 7H 1 OUTPUT / 1 HO . 9X . 2HP4 . 1 TX . 2HT4 . I 7X « 2HT 1 . 20X , 1 0A6/ 
X3CE17.8.2X )//7X«2HPl 

1 . 17X.2HUS. 17X.2HU2, 17X.2HP2//) 

IC=1 

159 IF(COMPUT.EQ.O) GO TO 160 
GO TO 165 

160 IF (REAL.EQ.O )GO TO 164 

161 CALL PERNS(P1 »T1 .GAMMA. MU. R.P2.U2.NUMP) 

GO TO 165 

164 CALL NORMAL (PI ,T1 ) 

165 CALL SOLUT (U3.P3.U2.P2. NUMT .NUMP » U « P ) 

IF(REAL.EQ.O) GO TO 170 

UP=SQRT(GAMMA*R*T1/MU)*(SQRT(1 •+ (P/Pl — 1 . ) # ( GAMMA+1 • )/<2.*GAMMA) ) ) 

UP=UP/30.48 

GO TO 175 

170 CALL FTLUP (U. UP.- 1 .NUMP .U2.UST0) 

175 CONTINUE 

I F (NS • EQ» 1 > WR I TE ( 6 . 4456 ) 

WRITE(6.207)P1 .UP.U.P 

207 FORMAT (4(E17»8.2X) ) 

IF(NS.EO.l) WR I TE ( 6 . 4456 ) 

4456 FORMAT!//) 

PLAST = P 1 

READ NEXT P 
READ (5. NORMS ) 

IFCP1 .EQ. PLAST) GO TO 1 
GO TO 159 
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END 
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SIBFTC PEREX DECK 

SUBROUT I NE PEREX ( P4 . U4 . GAMM A .R . T4 « MU « U3 . M , P3 ) 

DIMENSION U3 ( 30 ) . P3 ( 30 ) 

REAL MU 

GAMY=2.*GAMMA/ (GAMMA-1 . ) 

P4=P4*1 • 0 1 325E6 
DO 10 1 = 1. M 

P3 ( I ) =P4*< 1 .-(GAMMA-1 . ) /2 .* ( ( ( U3 ( I )-U4 )*30«48 J/SQRT ( GAMMA*R*T4/MU ) 
1 ) ) **GAMY 

10 P3( I )=P3( I>/1 .01325E6 
P4=P4/1 .01 325E6 
RETURN 
END 



SIBFTC ROGO LIST 

SUBROUTINE ROGO (PN , TM , NUMT , NCAPX . GAM, S0R4 ) 

C PERRY NEWMAN. EQUILIBRIUM THERMODYNAMIC PROPERTIES WITH DERIVATIVES 

80 FORMAT (2I5.2E14.8) 

81 FORMAT ( A6 .3I5.3E14.8) 

82 FORMAT (5E 14. 8) 

83 FORMAT (E14.8.2I5) 

84 FORMAT 1 4E1 4.8. 15 ) 

90 FORMAT <//4H MU=E1 5 . 8 , 2x ,3HP0=E 1 5.8 . 2X . 5HRH00=E1 5. 8 , 2X . 3HA0=E1 5.8 ) 
93 FORMAT ( / ( 8E16. 8 ) ) 

97 FORMAT ( // . 1 X6HA ( I , J ) . 2X2HI = 1 4 . 2X2H J= 1 4/ ( 8E1 6 • 8 ) ) 

98 FORMAT ( 1 X2HP=E 1 5» 8 « 2X28H1 00 I TERAT I ONS-NONCON VERGENT ) 

100 FORMAT ( 1 6H1 EXPANS I ON PHASE) 

101 FORMAT ( 75H0EQU I LI BR I UM THERMODYNAMIC PROPERTIES WITH DERIVATIVES 
1 IN REAL GAS SYSTEM//) 

999 FORMAT!/) 

1000 FORMAT ( 1 HI »57X2HT=I5, 1X1HK.//) 

1 001 FORMAT ( 6X3HLOG . 5X3HL0G . 7X1 HZ . 9X4HH/RT . 6X3HS/R . 6X3HL0G .6X6HDRH0DT . 5 
1 X6HDRH0DP , 3X4HCP/R . 6X4HCV/R . 5X5HGAMMA , 3X6HGAMMAE , 3X4HA/A0/5X4HP/P0 
2 . 3X8HRH0/RH00 . 34X2HNE . 5X7H ( T/RHO ) . 4X7H < P/RHO ) // ) 

1002 FORMAT (F9.1 ,F10.4,F9.4,F11.4,F10.4,F9.3.F10.3,F10.3,F10«3,F10.3,F8 
1 .3.F8.3.F9.3 ) 

REAL NO. M. LAMB, LAMBDA, MU. NE.NEGFRT.LOGNE.MASSFR.LM IN 

INTEGER F(30),V(30,10), PUN 

DIMENSION SPEC IE <30 ) » LB < 30 ) . M < 30 ) , DELHF ( 30 ) , BETA < 30) ,NDEBUG<30) 

1 . IPIVOT (10>»R(10.10> .SUMAY <10,1 ) , G < 30 . 30 ) . E < 30 ,30 ) , BE < 30 » 1 0 ) , ALPHA 
2E ( 30 ,10) .OMEGA <30. 10) .OMEGAX <30 , 1 0 ) » X0MEG<30.4) ,XOMEGX(30 

3.4 ) ,Z<30 ) .SIGMA < 1 0 > « U < 10 ) .DELTA < 10 ) .GAMMA < 1 0 ) » XX < 1 0 ) ,Q < 30 ) , Y < 30 ) , X 
4 (30) »A (30.9) . HORT ( 30 ) .FORT (30) » NEGFRT < 30 ) .PI <9,2 ) « XPRI ME <30 ) .MASSF 
5R ( 30 > ♦ CAPX ( 50 ) « Y I NT < 30 ♦ 1 ).CSUBP<30) .PS I (30,2) . CON < 1 0 < 2 ) , DXDT < 30 ) <R 
6R( 10. 10 > .0 ( 10. 1 0) , ABL < 30 ) * SBL ( 30 ) » HBL (30) ,PN<30) »TM(30) 

WR I TE (6 « 1 00 ) 

WRITE(6. 101 ) 

CALL SL I TET ( 1 , J J ) 

IF ( JJ.EQ.2 ) GO TO 7777 
H=6»625 1 7E-27 
XK=1 .38044E-1 6 
PREF =1.01 3250E + 6 
N0=6 .02322E+23 
C=2 .99793E+ 1 0 

1 READ (5 , 80 ) NUMSP, JINDX.EA.ER 
C IF NDEBUG EQUALS 0. DEBUG 

DO 3 1=1, NUMSP 

READ <5, 81 ) SPECIE! I ) ,LB < I ) ♦ F ( I ) , NDEBUG < I ) , M ( I ) , DELHF < I ) .BETA ( I ) 

I L =LB < I ) 
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READ (5, 82 ) (G(I,L),E(I.L).L=1.IL) 

I F ( F ( I ).EQ.O) GO TO 3 
IF(F(I ).EQ#2) GO TO 123 

READ (5.84 ) (BE ( I »L ) .AlPHAE ( I .L > .OMEGA ( I .L ) ♦ OMEGAX ( I *L ) » V( I »L ) »L = 1 . I 
1L > 

GO TO 3 

123 READC5.82) BE ( I . 1 ) , ALPHAE ( I . 1 ) 

READ (5. 82) <XOMEG( I .LW) .XOMEGX( I ,LW) ,LW=1 .4) 

3 CONTINUE 

READ (5. 82 ) MU , ( V I NT ( I . 1 ) . I = 1 .NUMSP) 

READ (5, 82) ( (A (I «J)«J=1 , JINDX) ,1=1 .NUMSP) 

7777 CONTINUE 

RH00=PREF*MU/(N0*X<*273. 15) 

AO=SQRT ( GAM * ( PREF/RHOO ) > 

WR I TE ( 6 , 90 ) MU.PREF.RHOO.AO 
DO 278 KP= 1 , NCAPX 

278 C APX ( KP ) = ALOG 1 0 (PN ( KP ) ) 

279 DO 19 <1=1 . NUMT 
T = TM ( K I > 

<T=T 

WRITE(6. 1000 > KT 
WR I TE (6 , 1001 ) 

NY = 1 

PART=H*C/ ( XK*T ) 

DO 9 1=1 .NUMSP 
I L=LB ( I ) 

I F (F ( I ).EQ.l ) GO TO 111 
I F ( F ( I ).EQ.2) GO TO 1 12 
QSUM=0. 

FPQSUM=0. 

SPQSUM=0. 

DO 2 L = 1 , IL 
2 ( L ) =PART*E ( I ,L) 

GEZ=G ( I , L ) *EXP ( -Z (L ) ) 

QSUM=QSUM+GEZ 
FPQSUM=FPQSUM+GEZ*Z ( L ) 

2 SPQSUM=SPQSUM+(Z(L)-2. )*GEZ*Z(L) 

FPQSUM=FPQSUM/T 

SPQSUM=SPQSUM/T**2 

Q I = (M ( I )*T*.32807618)**1 .5*QSUM*. 13623883*T 
GO TO 71 
111 QSUM=0. 

FPQSUM=0. 

SPQSUM=0. 

DO 11 L= 1 , I L 
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Z(L)=PART*E( I « L ) 

SIGMA(L)=PART*(BE< I * L >-.5*ALPHAE( I *L ) ) 

U(L)=PART* (OMEGA ( I • L ) -2 . *OMEGAX ( I « L ) ) 

DELTA (L)=ALPHAE { I «L > * 1 •/ (BE ( I « L ) - • 5* ALPHA E ( I « L ) ) 

GAMMA (L )= (BE ( I .L ) /OMEGA ( I ,L ) )**2*1 ./( 1 • -«5*ALPHAE ( I *L)/BE( I <L) • 

XX (L ) a OMEGAX ( I iL)/( OMEGA ( l i L ) -2 • *OMEGAX ( t<L)) 

THREE=0. 

FOUR=0. 

FIVE=0. 

NV=V ( I » L ) + 1 
DO 4 IV=1 »NV 
W= I V- 1 

CC=(1.-W*DELTA(L) ) 

A A=S I GMA ( L )*CC 

BB=U ( L )* ( W-XX (L >*W* (W-l . ) ) 

ONE1 *1 ./AA+8.*GAMMA(L >/ ( AA**2*CC >+ . 33333333+AA/l 2 . 

TW02= 1 ./AA+16.*GAMMA(L ) / ( AA**2*CC ) - AA/1 2 . -384 . *GAMMA (L ) **2/ ( AA**3* 
1CC**2 ) 

THREE = THREE+ ONE 1 *EXP ( -BB ) 

FOUR = FOUR+ (BB*ONE 1 +TW02 >*EXP ( -BB ) 

4 F I VE=F I VE+ ( (BB**2*0NE1+2.*B3*TW02+GAMMA (L)/ ( AA**2*CC ) * (48. -3456 . *G 
1 AMMA (L)/(AA*CC)+46080.*GAMMA (L ) **2/ ( AA**2*CC**2 ) 1+2./AA ) *EXP (-BB ) ) 
GEZ=G ( I .L)*EXP(-Z(L) ) 

Q ( L ) =THREE/BET A ( I ) *GEZ 
QSUM=QSUM+Q(L) 

FPQSUM=FPQSUM+(FOUR+THREE*Z(L ) ) *GEZ 
1 1 SPQSUM=SPQSUM+(Z(L)*(Z(L>-2. ) *THREE+2 . * ( Z ( L ) - 1 . ) *FOUR+F I VE ) *GEZ 
FPQSUM = FPQSUM/(T*BETA( I ) ) 

SPQSUM = SPQSUM/ ( T **2*BET A ( I ) ) 

Q ! - (M ( I >*T*. 328076 18 )**1 .5*QSUM*. 13623883*T 
GO TO 71 

112 S I GMA ( 1 ) =PART* ( BE (1*1 ) - • 5* ALPHAE ( I i 1 ) ) 

PROD= 1 • 

SUM 1 =0 • 

SUM2 = 0. 

DO 32 LW= 1 *4 

U(LW) =PART* (XOMEG ( I ,LW l-XOMEGX ( I *LW ) ) 

PROD=PROD* ( 1 • -EXP ( -U ( LW ) ) ) 

BTM=EXP(U(LW) )-l . 

SUM1=SUM1+U(LW)/BTM 

32 SUM2=SUM2+U(LW)**2*EXP(U(LW> )/BTM**2 
QSUM=G( I t 1 )/ (BETA ( I ) *SIGMA ( 1 ) *PROD ) 

FPQSUM= ( 1 ,+SUMl )*QSUM/T 

SPQSUM= (SUM1 **2+SUM2 )*QSUM/T**2 

Q 1 = (M ( I >*T*. 328076 18 )**1 .5*QSUM*. 1 3623883*T 
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CO 


71 HORT ( 1 >=2.5+T/QSUM*FPQSUM+DELHF ( I )/ <NO*XK*T ) 

FORT ( I >=DELHFC I > / ( NO*XK*T ) - ALOG ( Q I ) 

NEGFRT ( 1 )=-FORT ( I ) 

9 CSUBP ( I ) =2.5+2.*T/QSUM*FPQSUM- ( T*FPQSUM/QSUM ) **2+T**2*SPQSUM/QSUM 
NNN=1 

DO 39 KP= 1 » NCAPX 

9500 Pa ( 1 0»**CAPX (KP) >*PREF 
NYY=NY 

DO 500 I = 1 * N UM S P 
500 Y t I ) =Y I NT ( I *NY) 

9501 *NUMIT = 0***^* 

MM=UINDX+1 
DO 301 J=1 ♦ J I NDX 
0( J.MM )=0. 

DO 60 I=1«NUMSP 
60 0 ( J*MM )=0( J«MM )+A ( I « J )*Y< I ) 

0(MM« J)=0( J«MM) 

301 CONTINUE 

0 ( MM « MM ) =0 • 

50 CONTINUE 
YBAR=Y ( 1 ) 

NUM I T=NUM I T+ 1 

IFINUMIT.EQ. 101) GO TO 390 
DO 5 1=2* NUMSP 

5 YBAR=YBAR+Y ( I ) 

DO 7 K=1 ♦ JINDX 
R ( 1 «K)=0. 

DO 6 1 = 1 * NUMSP 

6 R ( 1 «K) =R ( 1 »K ) + A ( 1*1 )*A ( I *K )*Y ( I ) 

7 CONTINUE 
I COUNT = 1 
J J =2 

DO 28 J=JJ. JINDX 
DO 18 K=J, JINDX 
R(J«K)=0« 

00 8 1=1 « NUMSP 

8 R ( J.K ) =R ( J«K )+A (1 t J)*A( I *K ) *Y< I ) 

18 CONTINUE 

DO 1 0 K = 1 . I COUNT 

10 R( J.K)=R(K«J) 

1 COUNT = ICOUNT+ 1 
28 JJ=1+I COUNT 

DO 301 1 J=1 « MM 
R ( U « MM ) =0 ( J » MM ) 

3011 R ( MM « J ) =0 ( J « MM ) 
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DO 302 J=1 t MM 
DO 302 K=1 «MM 

302 RR(J.K )=R( J.K ) 

PYBAR=PREF*YBAR 
DO 304 J=1 .JINDX 
SUMAY ( J * 1 ) =0 • 

DO 303 I=1.NUMSP 
THIS=P*Y ( I J/PYBAR 
IFCTHIS.LE.O. > GO TO 303 

SUMAY ( J . 1 )=SUMAY( J. 1 )+A < I « J >*Y < I >* (FORT ( I )+ALOG(THIS ) ) 

303 CONTINUE 

304 CONTINUE 
SUMAY ( MM « 1 ) = 0. 

DO 305 1=1 .NUMSP 
TH I S=P* Y ( I l/PYBAR 
IFCTHIS.LE.O.) GO TO 305 

SUMAY (MM. 1 ) =SUM AY ( MM . 1 )+Y ( I >* (FORT ( I )+ALOG(THIS ) ) 

305 CONTINUE 
MN= 1 

NM AX= 1 0 

CALL SIMEO(R.MM. SUMAY. MN. DETERM. IPIVOT.NMAX. I SCALE ) 

DO 306 J=1 .JINDX 

306 PI ( J . 1 ) =SUMAY ( J ♦ 1 ) 

U=SUMAY ( MM . 1 ) 

LM I N= 1 . 

LCOUNT=0 

DO 40 1=1 .NUMSP 

AP I =0 • 

DO 401 J=1 .JINDX 

401 API=API+A ( I , J)*PI ( J. 1 ) 

TH I S=P* Y ( I ) /PYBAR 
IF(THIS.LE.l .E-38) GO TO 402 

X ( I ) =Y ( I ) * ( NEGFRT ( I > -ALOG ( TH I S ) +U+ 1 . + API ) 

GO TO 403 

402 X ( I ) =0 • 

403 I F ( X ( I ) ) 20.30.40 

20 LAMB =— Y ( I )/(X( 1 )-Y< I ) ) 

I F (LAMB .GT • 0 • ) GO TO 21 
Y ( I ) =0 • 

GO TO 50 

21 LCOUNT = 1 

LM I N=AM INI (LMIN.lAMB) 

GO TO 40 

30 IF (Y( I ) *EQ. 0 « >GO TO 40 
LCOUNT = 1 
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LAMB= 1 . 

LM I N=AM INI (LM1N.LAMB) 

GO TO 40 

40 CONTINUE 

I F (L COUNT .EQ » 0 ) GO TO 51 
LAMBDA= .999999*LMIN 
DO 41 1=1 .NUMSP 

41 Y ( I ) = C 1 .-LAMBDA )*Y ( I )+LAMBDA*X ( I ) 

GO TO 50 

51 DO 52 1=1. NUMSP 

I F ( Y ( I ) »EQ» 0 » ) GO TO 52 

IF ( ABS < X < I )-Y (I))/Y(I). GE.ER. OR. ABS C X ( I > — Y ( 1 ) > .GE.EA )G0 TO 53 

52 CONTINUE 
GO TO 29 

53 XBAR=0 . 

LAMBDA= 1 . 

LASTCT=0 

531 DO 54 1=1. NUMSP 

XPRIME ( I ) = ( 1 .-LAMBDA )*Y < I )+LAMBDA*X( 1 ) 

54 XB AR = XBAR+XPR I ME ( I ) 

DFLAMB=0. 

DO 541 1=1. NUMSP 

TH 1 S=P*XPR I ME ( I )/<PREF*XBAR) 

IF (THIS.LE.O. ) GO TO 541 

DFLAMB=DFLAMB+(X( I ) — Y ( I ) ) * ( FORT ( I >+ALOG(THIS) ) 

541 CONTINUE 

IF CDFLAMB.GT.O. )G0 TO 56 

542 DO 55 1=1, NUMSP 

55 Yd ) =XPR I ME < I ) 

GO TO 50 

56 LASTCT = LAST CT+ 1 

IF (LASTCT.EQ.4 )GO TO 542 
LAMBDA = . 9*LAMBDA 
XBAR=0 . 

GO TO 531 
29 XBAR=0. 

ONE=0 • 

TWO=0 • 

DO 57 1=1, NUMSP 
XBAR=XBAR+X< I ) 

ONE=ONE+X ( I ) *HORT ( I ) 

MASSFR ( I ) =X ( I )*M< I ) 

I F ( Y ( I J.EQ.O. ) GO TO 57 
TWO=TWO+X ( I ) * (FORT ( I )+ALOG(X( I) >) 

57 CONTINUE 


APPENDIX 



RECIPZ=1 ./(MU*XBAR) 

CAPU=CAPX <KP ) +ALOG 1 0 ( REC I PZ*273 . 1 5/T ) 

CHORT=MU*ONE 

SOR=CHORT -MU*(XBAR*ALOG(P/(PREF*XBAR) )+TWO) 

RHO=P/(XBAR*NO*XK*T ) 

NE=X( 1 )*RHO*NO 

IF(X(1 >.GT. (6.*10.**<-1 1 ) ) ) GO TO 58 
LOGNE=— 0 • 

GO TO 391 

58 LOGNE=ALOG 10 (NE ) 

59 GO TO 391 

390 WRITE (6 » 98 )P 
GO TO 39 

391 CONTINUE 

DO 600 I = 1 » NUMSP 
TEST=X( I )*P/(XBAR*PREF ) 

IFITEST.LT. 10. **(-20)1 GO TO 601 

PS I ( 1.1 )=X<I >/T*(HORT< I ) -FORT ( I ) -ALOG (TEST ) ) 

PSIII .2 ) =-X ( I ) 

GO TO 600 

601 PS I ( I . 1 ) =0 • 

PSI ( 1 .2 ) =0. 

600 CONTINUE 

DO 602 J=1 , JINDX 

CON ( J » 1 ) = A ( 1 ,J)*PSI (1,1) 

CON ( J , 2 ) =A ( 1 «J)*PSI (1.2) 

DO 603 1=2 .NUMSP 

CON( J* 1 ) =CON ( J » 1 )+A ( I * J)*PSI ( I • 1 ) 

603 CON ( J i 2 ) =CON ( J , 2 ) + A ( I »J)*PSt (1,2) 

602 CONTINUE 

CON (MM, 1 )=PSI (1,1) 

CON ( MM , 2 ) =PS 1 (1,2) 

DO 607 1=2, NUMSP 

CON (MM, 1 ) =CON ( MM , 1 )+PS! (1,1) 

607 CON ( MM , 2 ) =CON ( MM , 2 ) +PS 1(1,2) 

NC=2 

CALL SIMEQ(RR, MM, CON, NC, DETERM, IPIVOT,NMAX, ISCALE) 
DO 604 J= 1 , MM 
PI (J, 1 )=CON( J. 1 ) 

604 PI ( J , 2 ) =CON ( J , 2 ) 

DHDT=0. 

DO 605 1=1, NUMSP 
SUMAP=A ( I . 1 )*PI ( 1 , 1 ) 

DO 606 J=2, JINDX 
606 SUMAP = SUMAP+A ( I , J ) *P I ( J , 1 ) 


to 
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to 


DXDT ( I ) =PS I < I , 1 >-X ( I ) * < P I ( MM , 1 ) +SUMAP ) 

605 DHDT = DHDT+ { X ( I ) *C SUBP ( I > + T*H0RT ( I ) *DXDT ( I ) ) 

DRHODT = T*P I (MM* 1 )-l . 

DRH0DP= 1 • + P I ( MM ♦ 2 ) 

CPOR=MU*DHDT 

XZ=XBAR*MU 

CV0R=CP0R-DRH0DT**2/DRH0DP*XZ 
XGAMMA=CPOR/CVOR 
GAMMAE=XGAMMA/DRHODP 
RHOR= 1 ./ ( 1 0.**CAPU) 

TEST* (GAMMAE/GAM *P/PREF*RHOR > 

I F ( TEST »LE • 0 • ) WR I TE ( 6 ♦ 1 002 ) C APX ( KP ) , CAPU * XZ * CHORT « SOR *LOGNE«DRH 
1 ODT * DRHODP ♦ CPOR « CVOR ♦ XGAMMA , GAMMAE * AOAO 
AOAO=SQRT (TEST ) 

WRITE (6* 1002 ) CAPX(KP > , CAPU * XZ * CHORT * SOR *LOGNE . DRHODT * DRHODP , CPO 
1 R « CVOR .XGAMMA, GAMMAE*AOAO 
IF(NNN.EQ.IO) GO TO 608 
GO TO 609 

608 WRITE (6, 999) 

NNN = 0 

609 NNN=NNN+1 
ABL(KP)=AOAO*AO 
SBL(KP)=SOR 

HBL (KP ) =CH0RT*T*8 . 3 1 469E7/MU 
39 CONTINUE 

CALL INTER (PN ,T .NCAPX * SBL * HBL , ABL * S0R4 * K I ) 

19 CONTINUE 
RETURN 
END 
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SIBFTC INTER DECK 

SUBROUT I NE I NTER ( PN . TM , N « SBL . HBL . ABL . S0R4 * M ) 
DIMENSION PN ( 30 ) « SBL ( 30 ) *HBL ( 30 ) • ABL ( 30 ) • 

1 AM ( 30 ) « HM ( 30 ) 

COMMON/ BLOCK/HM . AM « M « P3 ( 30 ) 

DO 5 I * 1 * N 

5 PN ( I )=ALOG 1 0 (PN ( I ) ) 

NP«1 

CALL FTLUP < S0R4 «P3 ( M ) , NP , N. SBL ♦ PN ) 

NP=— 1 

CALL FTLUP t P3 ( M ) « HM ( M ) t NP * N « PN « HBL ) 

CALL FTLUP < P3 ( M ) » AM (M)t NP tNi PN < ABL ) 

DO 45 1=1 »N 
45 PN < I ) = 10«»*PN( I ) 

P3(M)=10t**P3(M) 

RETURN 

END 
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SIBFTC INTEG DECK 

SUBROUT I NE I NTEG ( DU , U4 > 

EXTERNAL FUNC 

DIMENSION AM ( 30 ) * HM ( 30 ) , DU ( 30 ) * SU ( 1 , 30 ) 

COMMON/BLOCK/HM i AMi MiP3 (30 ) 

DO 1 I =1 

1 AM ( I )-l #/AM( I ) 

DU ( 1 ) =U4 
DO 10 1=2, M 
NN= I 

CALL MGAUSS ( HM ( 1 ) , HM ( I ) , NN « SU (1,1) «FUNC,FOFx, 1 ) 
10 DU< I )=U4-SU( 1 , I )/30*48 
RETURN 
END 



SIBFTC PERNS DECK 

SUBROUTINE PERNS (Pi »T1 , GAMMA .MU .R .P2 * U2 . NUMP ) 

REAL MU 

DIMENSION P2 ( 30 ) * U2 ( 30 ) 

GAMM=GAMMA-1 . 

GAMP=GAMMA+ 1 . 

SQGRT = SQRT ( GAMMA*R*Tl /MU ) 

DO 10 1=1. NUMP 

U2 ( I ) =SQGRT # (P2 ( I l/Pl-l . >*SQRT ( ( 2 . /GAMMA 1 / < GAMP*P2 ( I )/Pl+GAMM) ) 
10 U2( I ) =U2 ( I 1/30.48 
RETURN 
END 


to 

UI 
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$ I BFTC NORMAL DECK 

SUBROUTINE NORMAL ( P 1 0 * T 1 0 ) 

C 

C P-886.5 

C NORMAL SHOCK PROGRAM 

C PROGRAMMED FOR THE IBM 7094 

C YIELDING SOLUTIONS FOR FLOW PARAMETERS IN ARBITRARY GAS 

C MIXTURES IN THE FOLLOWING SITUATIONS- 

C 1. BEHIND NORMAL SHOCK 

C 
C 

DIMENSION YSTO ( 30 ) . SHBL (8 ) 

C 

C EQUILIBRIUM INPUT 

REAL MU 

INTEGER PUN * COMPUT t REAL .EXPO 

COMMON /BLOCK 1 / I CODE ( 30 ) *F(30 ) .CAPM<30) . DHF 0 ( 30 ) «L ( 30 > ♦ 6 (30 *30 ) « 
1 SMLE ( 30 ♦ 30 ) * CAPLAM ( 30 « 30 ) * OMEG ( 5 . 30 • 30 > • A < 1 0 • 30 > • CONR . CONPRF ♦ 
2C0NN0 . CONH . CONK . P I .EPS1 .NIT.EPS2f I Cl 
COMMON PN ( 30 ) . TM ( 30 ) ♦ NUMT « NC APX . S0R4 *U4*P4.T4* U3 ( 30 ) » 

I PUN * EXPO, REAL, COMPUT, PI .T1 t P2 (30 ) ♦ NUMP . USTO ( 30 > tU2 (30 > 

COMMON I SPEC (30 ) . JMOL (10) • M * N * BET A ( 5 ) .NS* 

1 AMC ♦ NB « NBT A (5) . MU . GAMMA . R .NAME (10) 

C 

c 

C SHOCK PROGRAM INPUT 

C 
C 
C 

IF(NS.EQ.l) WR I TE ( 6 » 400 1 ) P10.T10 
4001 FORMAT ( 1 H 1 ♦ 1 3H NORMAL SH0CK/4H PI =E 1 5 • 7 . 2X ♦ 3HT 1 =E1 5 * 7 ) 

P 1 0 = P 1 0* 1 .01325E6 
DO 100 I J= 1 . NUMP 
CALL SL I TET ( 2 ♦ KK ) 

GO TO (44.45) , KK 

44 CALL TAPE(N. I SPEC * M . JMOL ) 

45 M = M 
US=USTO ( I J ) 

AMC=0. 

DO 52 I =1 .NB 
J=NBTA ( I ) 

52 AMC=AMC+BETA ( I ) *C APM ( J ) 

A 1 0=368 • *SQRT ( T 1 0/AMC ) 

AM10=US/Al 0 


C 


APPENDIX 



C SET UP GUESS Y< I) 

C 

IF < I SPEC (J ) «EQ • 1 )G0 TO 48 
DO 46 I =1 .N 
46 YSTO ( I ) = 1 .E-12 
GO TO 51 

48 YSTO ( 1 ) =0 • 

DO 49 1=2, N 

49 YSTOM >=YSTO(l 1+A fM, I 1*1 .E-12 
YSTO ( 1 ) =— YSTO ( 1 ) 

DO 50 1=2, N 

50 YSTO ( I )=1 .E-12 

51 CONTINUE 

DO 3035 1=1 ,NB 
J = NBTA ( I ) 

3035 YSTO(J)=BETA(I 1/AMC 

RHOl 0=P1Q/ (T 10/300. >*. 04061 9*AMC 

TP = T 1 0* ( 1 . + . 16* ( 1 . 5*AM 1 0**2+ 1 . ) /AM 1 0**2* ( AMI 0**2~1 . ) ) 

1 * ( 1 .-.5* (US-5000. 1/15000. 1 
RH02=RH010* (4.*AM 10**2 1 /< .5* AM10**2+2. ) 

C CONVERSION 

RHOl 0=RH01 0*1 • E-3 
RH02=RH02*1 .E-3 
US=US*30 • 48 

I F (T 1 0 .LE .800. ) CALL SLITE<4) 

CALL ECOM ( T 1 0 , P 1 0 , OOZ , HOZRT , HI 0 ,RHO 1 0 , YSTO , SOR 1 
C 

C STORE INITIAL P , T , RHO , AND U IN SHBL(l-4) 

C 

SHBL ( 1 >=P10 
SHBL ( 2 ) =T 1 0 
SHBL ( 3 ) =RHO 1 0 
SHBL ( 4 ) =US 

CALL SHOCK ( TP, SHBL, H 10, YSTO, RH02, US 1 
C 

C UPON RETURN SHBL ( 5-8 ) CONTAINS P2 , T2 , RH02 , AND UF 

C 

PPR INT=SHBL ( 5 ) /CONPRF 
C 

UPRINT=SHBL (81/30.48 
C 

P2 ( I J 1 =PPR I NT 
U2( I J 1 =UPR I NT 
IF(NS.NE.l 1 GO TO 100 

WRITE (6, 202 1 US , PPR I NT , SHBL ( 7 1 , OOZ , HOZRT , SOR , SHBL ( 6 1 , AMC 


i 
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202 FORMAT <//4H US=E1 5. 7//9X . 1 HP » 1 3X * 3HRH0 * 1 4X * 3H 1 /2 * 1 4X *5HH/ZRT . 
X12XOHS/R. 

1 14X* 1HT* 16X«2HM1//(7E17.8) ) 

WRITE <6*221 ) 

221 FORMAT ( //23H FINAL Y FROM I TERAT I ON . 4X ♦ 7HSPEC 1 ES// ) 

WRITE (6 *220) <YSTO< I ) * I CODE < I ) . 1 = 1 *N) 

220 FORMAT <E17*8«10X. 1 A6 ) 

100 CONTINUE 

P10=P10/1 .01325E6 

RETURN 

END 
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* IBFTC TAPE DECK 

SUBROUTINE TAPE (N , I SPEC . J , JMOL ) 

C 

COMMON /BLOCK 1 / 1 CODE ( 30 ) »F (30 ) *CAPM (30 ) *DHFO ( 30 ) *L ( 30 ) *G <30*30 ) * 
1 SMLE ( 30 * 30 ) * CAPLAM ( 30 * 30 ) * OMEG (5*30*30 ) * A ( 1 0*30 ) «CONR*CONPRF • 
2C0NN0 * CONH »CONK *P I ,EPS1 *NIT*EPS2* IC1 
D I MENS I ON BLOCK (150)* LBLOCK ( 30 ) . I SPEC ( 30 ) * UMOL (10)* 


1 OBL ( 5 * 30 ) 08967 

READ (9) (LBLOCK ( I) . I =1 .30 > 08967 

DO 1 I C = 1 . N 08967 

I SP = I SPEC ( I C ) 08967 

1 I CODE ( I C ) =LBLOCK (ISP) 08967 

C 08967 

READ (9) (Block (I >. 1=1 *30> 08967 

DO 2 IC = 1 *N 08967 

I SP= I SPEC ( I C ) 08967 

2 F< IC)=BLOCK( ISP) 08967 

C 08967 

READ (9) (BLOCK ( I ). 1=1 *30) 08967 

DO 3 IC= 1 * N 08967 

I SP= I SPEC ( I C ) 08967 

3 CAPM { I C ) “BLOCK (ISP) 08967 

C 08967 

READ (9 ) (BLOCK ( I ) . 1 = 1 ,30 ) 

DO 4 IC=1 *N 08967 

I SP= I SPEC ( I C ) 08967 

4 DHFO ( I C ) “BLOCK (ISP) 08967 

C 08967 

READ (9) ( LBLOCK ( I) « 1=1 *30) 08967 

DO 5 I C = 1 * N 08967 

I SP= I SPEC ( I C ) 08967 

5 L( IC) “LBLOCK (ISP) 08967 

C 08967 

IC = 1 

DO 6 1=1*30 08967 

READ (9) (BLOCK (IL) . IL=1 *30) 

IF ( 1 SPEC ( I C ) — I )6*55*6 

55 DO 56 L I =1*30 08967 

56 G ( L I • I C ) =BLOCK ( L I ) 

IC=IC+1 

6 CONTINUE 08967 

C 08967 

IC = 1 

DO 7 1=1,30 08967 

M READ (9) (BLOCK(IL) . IL=1 *30) 08967 

co 
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08967 


1 F ( I SPEC ( I C ) -I ) 7 , 65 , 7 

65 DO 66 L I =1 .30 

66 SMLE (LI , I C ) = BLOCK (LI) 

IC=IC+1 

7 CONTINUE 08967 

IC=1 

DO 12 1=1.30 

READ (9) (BLOCK (IL) » IL=1 .30) 08967 

I F ( I SPEC (IC)-I >12.13.12 

13 DO 125 LI =1.30 08967 

125 C A PL AM ( L I , I C > =BL0CK (LI) 

IC=IC+1 

12 CONTINUE 08967 

C 08967 

I I C= 1 

DO 8 1=1 .30 08967 

READ (9 > ( ( OBL ( I C . I L > . I C= 1 . 5 > . IL= 1 . 30 > 

I F ( I SPEC ( IlO- I >8,75.8 

75 DO 76 L 1=1 .30 

DO 76 I C= 1,5 08967 

76 0MEG< IC.LI , I IC)=0BL( IC.LI > 

I I C= I IC + 1 

8 CONTINUE 

C 08967 

IC=I 

DO 10 1=1.30 08967 

READ (9) ( BLOCK ( I J> , I J=1 , 10) 08967 

I F ( I SPEC (10-1)10.85,10 

85 DO 9 IJ=1,J 08967 

IJM=UMOL(IJ> 08967 

9 A( IJ.IC) =BLOCK ( I JM ) 

IC= IC+1 

10 CONTINUE 08967 

C 08967 


C0NR=8 .31 46938E7 
CONPRF = 1 .01325E6 
C0NN0=6 .02322E23 
C0NH = 6 » 625 1 7E— 27 
CONK= 1 . 38044E- 1 6 

PI=3. 14159 08967 

NIT = 300 
EPS1 = 1 » E— 7 
EPS2= » 0 1 

I F ( I SPEC ( N ) -28 114,15,14 
14 10=1 
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SIBFTC SHOCK DECK 

SUBROUT I NE SHOCK ( TGUESS . BLOCK , H l 0 . YSTO » RH02 * US ) 

C 

C THIS SUBROUTINE USES A ONE-DIMENSIONAL NEWTON-RAPHSON ITERATION 

C SCHEME TO FIND TEMPERATURE AND PRESSURE AT EQUILIBRIUM BEHIND 

C INCIDENT SHOCK. IT WILL CALL SUBROUTINE ECOM TO COMPUTE THE 

C EQUILIBRIUM PROPERTIES. 

DIMENSION YSTO (30) 

DIMENSION T (2 ) .H (2 ) .BLOCK (8 ) 

C 

C 

C LET RH02 = FIRST RHO 

C 

VEL1 (AA)=C*D/AA 

PRES 1 (AA«BB)=B+C*D**2-AA*BB**2 
ENTH1 ( A A )=H1 0+ (D**2 )/2.-< AA**2 )/2. 

DELT =10* 

I T = 3 

EPS5=1 .E-5 
NCOUNT= 1 
B=BLOCK ( 1 ) 

C=BLOCK (3 ) 

DEBLOCK (4 > 

U2=VEL1 (RH02 ) 

P2=PRES 1 (RH02.U2) 

H2=ENTH 1 (U2 ) 

C 

C COMPUTE FIRST POINT 

C 

I TT = I T 

7 T ( 1 )=TGUESS 

CALL ECOM (T ( 1 ) . P2 «OOZ . HOZRT * H ( 1 ) .RHO. YSTO .SOR ) 

C 

C COMPUTE SECOND POINT 

C 

T ( 2 ) =T ( I ) +DELT 

CALL ECOM (T (2 ) . P2 . OOZ . HOZRT . H ( 2 ) .RHO ♦ YSTO ♦ SOR ) 

S= ( H ( 2 ) — H ( 1 > )/ ( T ( 2 ) -T ( 1 ) ) 

T ( 1 ) =T ( 2 ) 

C 

C TEMPERATURE FROM FIRST ITERATION 

C 

T (2)=T (2)+ (H2— H (2 ) )/S 
HC 1 )=H(2 ) 

1 F (T (2 ) >25.25.8 


APPENDIX 


8 CALL ECOM(T<2> . P2 .OOZ .HOZRT ♦ H ( 2 ) .RHO « YSTO .SOR ) 

C 

C S IS SLOPE (H2— HI ) / ( T2-T 1 ) 

C 

85 S=(H(2)-H(1 ) )/(T(2)-T(l )) 

C T3 

T ( 1 ) =T ( 2 ) 

C 

c temperature from second ITERATION 
c 

T ( 2 ) =T ( 2 ) + ( H2-H < 2 > )/S 
H ( 1 ) =H ( 2 ) 

I F ( T ( 2 I >25.25.9 

9 CALL ECOM (T (2 ) .P2.00Z.H0ZRT .H (2 ) «RHO. YSTO.SOR ) 

C 

C IF ITT IS GREATER THAN 2. ITERATE AGAIN ON TEMPERATURE WITH 

C FIRST PRESSURE 

C 

I F ( ITT-2U0.10.il 

10 SLAST=(H(2)-H(1 ) >/<T(2)-T(l > ) 

TLAST = T (2 > 

GO TO 12 
C 

11 I TT = I TT— 1 
GO TO 85 

C 

C 

C TEST RHO FOR CONVERGENCE 

C 

12 I F ( ABS ( ( RH0-RH02 ) /RH02 ) — EPS5 >20.20. 13 
C 

C NON-CONVERGENCE- 

C 

C COMPUTE NEW PRESSURE AND CONTINUE ITERATION ON TEMPERATURE AND 

C PRESSURE UNTIL RHO CONVERGES 

C 

13 RH02 =RHO 
U2=VEL 1 ( RH02 > 

P2=PRES 1 CRH02.U2) 

H2=ENTH 1 ( U2 ) 

NCOUNT = NC0UNT + 1 

T ( 1 ) = TLAST 

14 CALL ECOM (T ( 1 ) , P2 . OOZ . HOZRT . H < 1 ) .RHO. YSTO.SOR) 

145 S=SLAST 

15 T(2)=T( 1 )+(H2-H(l> )/S 
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CALL ECOM (T (2) .P2 * OOZ « HOZRT . H (2 ) *RHO « YSTO .SOR ) 
SLAST = ( H ( 2 > -H ( 1 ) )/(T(2 )-T(l ) ) 

TLAST = T (2 ) 

GO TO 12 
C 

C CONVERGENCE - STORE OUTPUT 

C 

20 U2=VEL1 (RHO ) 

UF=US— U2 
RHO=RHO*l «E3 
BLOCK (5 )=P2 
BLOCK ( 6 ) =TLAST 
BLOCK ( 7 ) =RHO 
BLOCK (8)=UF 
RETURN 
C 

C TEMPERATURE ESTIMATE TOO HIGH - ADJUST 

C 

25 TGUESS= (TGUESS—Tl 0 )/2» 

GO TO 7 
ENO 


to 

co 
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GO 

SIBFTC ECOM DECK 

SUBROUT I NE ECOM ( T * PSTO , OOZ , HOZRT , H , RHO ♦ YSTO , SOR ) 

C 

C SUBROUTINE WHICH, GIVEN A TEMPERATURE AND PRESSURE, COMPUTES 

C THE THERMODYNAMIC EQUILIBRIUM PROPERTIES OF A GAS DESCRIBED BY 

C THE INPUT. 

C 

REAL mu 

INTEGER PUN, COMPUT, REAL, EXPO 
DIMENSION SMALE (30 , 30 ) » X < 30 ) ,YSTO(30 ) 

D I MENS ION E ( 30 ) . Y ( 30 5 , Q ( 30 ) . CAPF I(30)»R<10.10>.B<10). 

3TEMPS (10) , BSUM ( 1 1,1) , ABLOCK ( 11,11) »PTEMP(30 ) «ZETA(30) . 

4ZETAPR (30 ) , ALAM (30 ) * 

5 IP I VOT (11 ) , DQ I NT ( 30 ) , Q I NT (30,30) 

C 

COMMON /BLOCK 1 / 1 CODE ( 30 ) , F ( 30 ) , CAPM <30 ) , DHFO ( 30 ) ,L ( 30 ) .G ( 30 * 30 ) » 
1 SMLE ( 30 , 30 ) , CAPLAM (30.30 ) , OMEG ( 5 , 30 . 30 > , A ( 1 0 , 30 ) ,CONR .CONPRF . 
2C0NN0 , CONH , CONK ,P I ,EPS1 .NIT.EPS2, 1C1 
COMMON PN ( 30 ) , TM ( 30 ) , NUMT , NCAPX , SOR4 , U4 , P4 , T4 *U3 (30) • 

1 PUN, EXPO, REAL, COMPUT, Pi *T1 ,P2 (30 > « NUMP , USTO ( 30 > »U2 (30 ) 

COMMON I SPEC (30 ) , JMOL ( 1 0 ) , M , N , BETA ( 5 ) ,NS« 

1 AMC»NB* NET A ( 5 ) , MU, GAMMA , R , NAME (10) 

C 

EQUIVALENCE (SMLE ( 1,1) , SMALE ( 1 , 1 ) ) , ( I CODE ( 1 ) ♦ CODE (1 ) ) 

C 

PI=3. 14159 
C = 2 » 9979 3E 1 0 
NCOUNT = 0 
LTEST=LTEST 
N2 = N 

DO 5 I = 1 ,N 
5 Y( 1 ) = YSTO ( I ) 

P=PSTO 

34 TK=CONK*T 
RT = CONR*T 

346 YBAR=0 • 0 

DO 347 1=1 .N 

347 YBAR=YBAR+Y( I ) 

DO 40 1=1 ,N 
TEMPI =0 
LEND=L ( I ) 

DO 37 L 1 = 1 .LEND 
I F t F ( I ) ) 31 .35.31 
31 PR0D=1 • 

DO 33 IC=1.1C1 
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c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 

c 


IF(OMEG< ICtLl * I )) 32 03*32 

32 PR0D=PR0D*< 1 .-EXP ( -C0NH*C*0MEG ( I C «L 1 « I )/TK) ) 

33 CONTINUE 

PART = ( T/ ( C A PL AM (LI* I )*PROD) >**F(I ) 

GO TO 36 

35 PART = 1 « 

36 QINT CL1 * I )=PART*G(L1 , I ) *EXP (-CONH*C*SMALE (LI t I >/TK> 

37 TEMPI =TEMP1+QINT(L1 * I ) 

Q ( I ) = (SQRT(2.*PI/C0NH*TK/<C0NH*C0NN0)*CAPM< I ) )**3>*TK/C0NPRF*TEMPl 
IF (Y( I ) ) 38 * 38* 39 

38 CAPFI ( I ) =0 
GO TO 40 

39 CAPFI ( I ) = Y < I ) * ( ALOG ( P/CONPRF ) + ALOG ( Y ( I ) /YBAR ) -ALOG ( Q ( I ) )+DHFO< I ) 
1/RT) 

40 CONTINUE 

CALL SLITET(4*JJ) 

GO TO (95*396) * J J 
396 DO 50 J=1 « M 
DO 50 K=1 *M 
R <K* J )=0*0 
B ( J ) =0 • 0 
DO 50 1 = 1 * N 
B ( J ) =B ( J ) +A ( J * I )*Y( I ) 

50 R(K* J ) =R (K* J )+A ( J, I )*A (K* I )*Y ( I ) 

SET UP MATRIX FOR SOLUTION OF EQUATIONS 

DO 60 J=1 *M 
TEMPS ( J ) =0 • 0 
DO 55 1=1 *N 

55 TEMPS ( J )= TEMPS ( J ) + A ( J * 1 )*CAPFI ( I ) 

BSUM ( J * 1 ) =B ( J ) + TEMPS ( J ) 

CONSTANT TERMS IN BSUM BLOCK 

DO 56 K = 1 *M 
K 1 = K + 1 

56 ABLOCK( J*K1 )=R(K* J) 

PI TERMS IN ABLOCK IN COLUMNS 2 THROUGH N+l 
60 ABLOCK ( J * 1 )=B( J) 

(X/Y) TERMS IN FIRST COLUMN 


CO 

CJ1 
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Ml =M+1 

ABLOCK ( Ml .1 )=0.0 
DO 61 K = 1 .Ml 
K 1 =K+ 1 

61 ABLOCK (Ml ,K1 )=B (< ) 

BSUM (Ml .1 1=0.0 

DO 62 I =1 ,N 

62 BSUM ( M 1 ,1 1 =BSUM (M 1 « 1 1+CAPFI (I ) 

C 

C MAT1NV EXPECTS AN M+l BY M+l MATRIX 

C 

CALL S I MEQ( ABLOCK ( 1 . 1 ) .Ml .BSUM (1.11.1 .DETERM. IPIVOT. 11.01 
C 

C RETURN WITH ANSWERS IN BSUM 

C 

ZETAP=BSUM (1.1 ) *YBAR 
ZERO=0. 

NEG=0.0 
DO 70 1=1 .N 
PTEMP ( I 1=0.0 
DO 65 J=1 .M 
J1=J+1 

65 PTEMP ( I ) =PTEMP ( I 1+BSUM( Jl. 1 1 *A < J ♦ I 1 *Y ( I ) 

ZETA ( I 1=-CAPF1 ( I 1+Y( I }*BSUM (1,1 l+PTEMP ( I 1 

C 

C TEST FOR NEGATIVE OR ZERO ZETA 

C 

68 IF (ZETA( I 1 169,695.70 

69 P I ECE=- Y ( I >/(ZETA(I 1— Y( I 1 1 
IF (PIECE1691 ,692,691 

691 NEG=NEG+ 1 
ALAM(NEG)=PIECE 
GO TO 70 

692 Y( I }=0 
ZERO= 1 . 

GO TO 70 

695 1F(Y( I 1 169.70.69 

70 CONTINUE 
C 

C FIND GREATEST NEGATIVE ZETA-Y 

C 

IF (ZERO 1700, 700. 698 

698 I F ( NCOUNT-N IT} 699 ,100,1 00 

699 NCOUNT = NCOUNT + 1 
GO TO 346 
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00 

-a 


700 I F ( NEG— 1 ) 78 ,71, 73 

71 ALAMPR=.999999»ALAM ( 1 ) 

GO TO 745 

73 ARG1 =A|_AM ( 1 ) 

DO 74 I =2, NEG 

72 ARG2=ALAM(I) 

ARG1 =AM INI (ARG1 , ARG2 ) 

74 CONTINUE 
ALAMPR=.999999*ARG1 

745 IIC=0 

75 ZETAP=0 

DO 76 1=1 ,N 

ZETAPR ( I ) = Y ( I ) +ALAMPR* ( ZETA ( I) -Y ( I) ) 

76 ZETAP=ZETAP+ZETAPR( I ) 

DLAM=0 

DO 77 1 = 1 ,N 
I F (ZETAPR ( I 1 177.77,765 

765 DLAM = DLAM+ (ZETA (1 ) — Y < I ) ) * ( ALOG (P/CONPRF l-ALOG ( Q ( I ) )+DHFO(I l/RT+ALO 
1 G ( ZETAPR < I 1 /ZETA P 1 1 

77 CONTINUE 
IFfDLAM >81 ,81,80 

80 I F ( I I C-3 ) 805 ,81 ,81 
805 IIC=IIC+1 

AL AMPR= ALAMPR* . 9 
GO TO 75 

78 A|_AMPR= 1 • 

GO TO 745 

C 

C CONVERGENCE TEST FOR Y(I)S 

C 

81 IF(ALAMPR-. 50)83, 815, 815 

815 DO 82 I = 1 , N 

IF (ZETAPR ( 1) 1813.816,81 3 
813 REL=Y ( I ) — ZETAPR ( I 1 

IF(ABS(REL)-EPS1 1818,818,83 
818 REL = ZET APR (I )/Y ( I 1- 1 • 

IF ( ABS ( REL ) -EPS 2 1 82 , 82 , 83 

816 I F ( Y ( I ) 1817,82,817 

817 GO TO 83 

82 CONTINUE 
C 

C Y ( I IS CONVERGE 

C 

DO 800 1=1 ,N 
800 Y ( I )=ZETAPR ( I 1 
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C 

C 

c 


GO TO 95 


NON-CONVERGENCE OF Y(I)S 


83 NCOUNT =NCOUNT+ 1 

I F ( NCOUNT— N IT) 84 . 1 00 « 1 00 

84 DO 85 1 = 1 »N 

85 Y( I )=ZETAPR( I ) 

C 

C REPEAT WITH NEW Y(I)S AND NO. OF ITERATIONS LESS THAN NIT 

C 

GO TO 346 
95 DO 201 1 = 1 .N 

201 X( I >=Y( I )*CAPM( I ) 

YBAR=0 .0 
CAPM 1=0 
DO 2026 1=1 .N 
YBAR=YBAR+Y ( I ) 

2026 CAPM I =CAPM I +X ( I )/CAPM( I ) 

CAPM 1=1 • 0/CAPM I 
Z=AMC/CAPMI 

ESUM=0 

DO 2029 1 = 1 .N 
QSUM=0 
DQ I NT ( I )=0 
LEND=L ( I ) 

DO 2028 LI =1. LEND 
SUM = 0 

DO 2027 IC=1 , IC1 
HOOTK=CONH*C*OMEG( I C*L 1 « I )/TK 
I F ( OMEG ( I C * L 1 . I ) >2000.2027.2000 
2000 SUM=SUM+HOOTK/ ( EXP ( HOOTK ) — 1 . ) 

2027 CONTINUE 

DQINT ( I ) =DQI NT ( I ) +QI NT (L 1 . I ) * ( F ( I ) /T* ( 1 .+SUM ) +SMALE ( L 1 . I )*C0NH*C 
1/<TK*T) ) 

2028 QSUM=QSUM+Q I NT <L1 . I ) 

E ( I ) = 1 ./CAPM < I )*< 1 .5*RT+RT*T/QSUM*DQINT ( I ) + DHF0 < I ) ) 

2029 ESUM=ESUM+X ( I ) *E ( I ) 

H0ZRT=CAPMI*ESUM/RT+1 .0 

H=HOZRT*CONR*T*Z/AMC 

TK=T*CONK 

FSUM=0 

DO 2040 I = 1 « N 

2033 IF(Y( I ) >2034.2034.2035 

2034 CAPFI ( I ) =0 
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GO TO 2040 

2035 CAPFI ( I ) =Y ( I >* ( ALOG (P/CONPRF ) +ALOG < Y < I ) /YBAR ) -ALOG (Q < I > >+DHF0< I > 
1/RT) 

2040 FSUM=FSUM+CAPFI ( I ) 

SOZR=HOZRT-CAPM I*FSUM 
SOR=SOZR*Z 

rho=p*capmi/rt 

U=CAPX+.43429*AL0G(273. 16/(Z*T ) ) 

OOZ=l ,0/Z 
DO 300 1 = 1 *N 
YSTO ( I ) =Y ( I ) 

300 X< I )=X( I)*CAPMI/CAPM( I ) 

RETURN 

100 WR 1 TE (6 *5000 > 

5000 FORMAT! 1 HO* 25H THIS CASE NON-CONVERGENT ) 

CALL EXIT 
END 


I 


co 
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SIBFTC FOFX DECK 

SUBROUTINE FUNC < DUM • FOFX > 
COMMON/BLOCK/HM « AM • M *P3<30 ) 
DIMENSION HM (30 ) * AM ( 30 ) 

CALL FTLUP ( DUM * A A « -2 ♦ M ♦ HM ♦ AM ) 

fofx=aa 

RETURN 

END 
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SIBFTC SOLUT DECK 

SUBROUTINE SOLUT (U3 *P3. U2 . P2 . M ♦ N .UR . P ) 
DIMENSION U3<30).P3C30) .U2 < 30 ) .P2 ( 30 ) .U < 2 ) 
FUNCD (PP.UU.R) =PP-UU*R 
FUNAB (P .PP.U.UU > = (P-PP)/(U-UU) 

C 

C USE END POINTS FOR FIRST INTERSECTION 

C 

MR = 1 
NR= 1 

1 IF (P2 (1 ) .GT.P2 t 2 > > NR=— NR 

! IF(P3(1 >.GT.P3(2> > MR=-MR 

P3I =P3 ( 1 ) 

P32=P3 t M ) 

P21 =P2( 1 ) 

P22=P2(N) 

U21=U2( 1 ) 

U22=U2(N> 

U31 *U3 ( I > 

U32=U3(M > 

5 AA=FUNAB ( P22 .P2 1 * U22 . U2 1 ) 

BB=FUNAB (P32.P31 .U32.U31 ) 

' CC=FUNCD(P2l .U21 ♦ AA ) 

DD=FUNCD(P31 .U31 *BB) 

UR= ( CC-DD > / t BB-AA ) 

PR=CC+UR*AA 

CALL FTLUP ( PR ♦ U ( 1 ) * NR . N *P2 ♦ U2 ) 

CALL FTLUP (PR * U ( 2 ) .MR.M.P3.U3) 

I F ( ABS ( (U ( 1 )-U<2) )/U(l > >-.0001 >12.12*10 
10 P31=P32 
P32=PR 
P21=P22 
P22=PR 
U31=U32 
U32=U ( 2 ) 

U2 1 =U22 
U22=U( 1 ) 

GO TO 5 
12 P=PR 

! RETURN 

END 

I 

I 
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Pressure, p^/p 


Figure 10.- Test gas shock speed as a function of buffer gas pressure for test gas of 90 percent N 2 and 10 percent CO 2 
and helium driver conditions of pWpg = 100 and T 4 = 300° K. 
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Figure 11.- Test gas shock speed as a function of buffer gas pressure for test gas of 9( 
hydrogen driver conditions of pWp n = 100 and Li = 3( 


